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ABSTRACT 


A three-dimensional (3D) coupled normal mode model for studying 
sound propagation in a complex coastal environment is developed. 
This development corresponds to a significant upgrade of an earlier 
version of the model in which a flat, rigid bottom was used. By 
imposing the general boundary conditions for an irregular, non- 
rigid bottom, the coupling coefficient integrals in the system of 
differential equations governing the mode amplitude are re- 
formulated. The model upgrade entails a numerical implementation of 
the revised formulae. With the improved physics, this latest 
version iS capable of modeling the 3D acoustic wave-field in 
shallow water where sound speed, water depth and sediment 
properties can vary with horizontal location. To demonstrate this 
enhanced capability, the model is used here to simulate the 
interactions of the normal modes as they propagate up a sloping 
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i. INTRODUCTION 


A. BACKGROUND 

There are three approaches to model three-dimensional (3D) 
sound propagation in the ocean: ray theory, parabolic 
equation approximation and normal mode theory. 

Ray theory gives an approximate, planewave-like solution 
to the wave equation, which is valid at high enough 
frequencies and in media with gradual variations. The ray 
solution is constructed by raytracing. The acoustic rays 
provide for a visual, physical description of sound 
transmission in the ocean. The ray solution, however, 
neglects sound diffraction and dispersion and thus needs 
corrections near caustics and CUtmiiges points. These 
corrections may sometimes be mathematically complicated. The 
Hamiltonian Acoustic Ray Tracing Program for the Ocean (HARPO) 
is the only 3D ray theory model available today. This computer 
code was originally developed by Jones et al. [Ref. 1] for the 
Gemputation of 3D rays. 

The parabolic equation approximation method (PE) was 
introduced by Tappert [Ref. 2]. PE is a "full-wave" method 
piteaccOUnts £Lor both sound diffraction and dispersion. It 
provides for numerical solutions to the wave equation which 


are accurate for energy propagating at low grazing angles. The 


accuracy generally degrades as the angle increases. The 
backscattered energy 1s generally neglected in this 
approximation. A versatile 3D PE model has been developed by 
Lee et al. [Ref. 3] using an implicit finite difference 
scheme. Another 3D PE model was developed earlier by Baer 
[Ref. 4] which uses a split-step Fourier algorithm. The PE 
model of Lee has a wider angle capability, i.e., it models 
sound energy travelling at steeper angles more accurately. 
Finally, normal mode theory describes sound propagation as 
a collection of eigenfunctions, called normal modes, which are 
a natural set of vertical vibration modes in the sound 
channel. Just like PE, normal mode theory is a "full-wave" 
approach. The original normal mode theory assumes a 
horizontally stratified propagation medium. This assumption is 
valid for many short-range, deep-water cases, where range and 
azimuthal variations are negligibly small. Pierce [Ref. 5] 
extended the theory to account for horizontal sound speed, 
bathymetry and bottom-property variations. These variations 
produce mode coupling phenomena (in which energy exchange 
between modes takes place). A 3D coupled normal mode model has 
been developed by Chiu and Ehret [Ref. 6]. This model is 
capable of simulating mode-mode interactions due to a 3D 
varying sound speed field. The effects of bottom bathymetry 


variations and sediment property, however, are not modeled. 


B. THESIS OBJECTIVES AND OUTLINE 

The main objective of this thesis is to improve the Chiu- 
Ehret [Ref. 6] 3D coupled normal mode model by including the 
effects of bathymetry variations and sediment properties on 
sound propagation. The improved model is useful for studying 
sound propagation in shallow water environments where 
Significant bottom interaction is expected. 

In Chapter II, 3D coupled mode theory is first reviewed. 
The formulae for the mode coupling coefficients in the system 
of differential equations governing the mode amplitude 
functions are derived. In the derivation, the general boundary 
conditions for an irregular, non-rigid bottom are used. 

In Chapter III, alternative expressions for the mode 
coupling coefficients are derived. These expressions allow for 
an easier numerical implementation. The improved model is used 
to examine the effects of a sloping bottom on upslope sound 
propagation. The validity of the adiabatic approximation is 
also examined. Conclusions are given in Chapter IV. 

The products coming out of this thesis are computer 
Subroutines to include bathymetry variations and sediment 
properties in the 3D coupled mode model of Chiu and Ehret 


[Ref. 6]. The new routines are listed in the Appendix. 


oe es 3D COUPLED NORMAL MODE THEORY 

In the mathematical Formulation that Follows, a 
cylindrical coordinate system will be used (see Fig. 1). The 
z-aX1iS 1S perpendicular to the ocean surface and iS positive 
downward, r is range from the source location (i.e., the 
origin) and @ is the azimuthal angle (positive clockwise). 
Sound speed in the water column, c,, is a function of r, z and 
6, where sound speed in the sediment, c,, is assumed to be r 
and @ dependent only. The density of the water column, p,, is 
considered to be constant. The density of the sediment, p,, is 
also considered to be constant. The water-sediment interface 


is located at z=H(r,@). 


A. THE MATHEMATICAL PROBLEM 
In the case of isodensity layers, the 3D, homogeneous, 
monofrequency Helmholtz Equation governing the acoustic 


pressure, p, is: 


V2piize, Oita 0) re) an) 


where k(z,r,9)=w/c(z,r,6@) is the acoustic wavenumber, w is the 
source angular frequency and c is sound speed (c, in the water 
layer and c, in the sediment layer). 

A quasi-separable solution to Eq. (1) is postulated, which 


is locally a linear combination of normal modes or depth 
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Figure 1. The model geometry 
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where R, are the mode amplitude functions and n is the mode 


number. 


The normal modes Z, are required to satisfy the depth 


equation: 





+ R°(2)r,0) = k* (ane) |Z 402, 2 0) eee (3) 
OZ 


where k, is the horizontal component of the wavenumber 


n 
(eigenvalue) associated with the n"” mode. 
It can be easily shown, using the boundary conditions for 


Z, {to be derived next) and the depth equation (Eq. (3)), that 


Nn 
the normal modes form a complete set of orthogonal functions, 
with the inverse of the medium density p as a weighting 


fFUNCENON In the normal? zace von. 
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where 6,, 1S the Kronecker delta. Note that the integration is 


Carried over the entire depth from 0 to ». Also, note that the 


density p in the model is considered to be constant in each 


layer (p, in the water and p, in the sediment). 


B. THE BOUNDARY CONDITIONS 
The appropriate boundary condition for the acoustic 


pressure at the sea surface is 


Pp, t2=077 0). = 70 (5) 


where the subscript 1 denotes the water column. This pressure 
release condition implies that the normal modes Z, must also 


be zero at the sea surface, l.e., 


Zo 7Z=07,0) = 6 (6) 


It also implies that the horizontal derivatives of Z, at z=0 


are zero, i1.e., 


d2( 2-08, 0) A (7) 
or 
1 02Z,(2=0;7r,9) =, (8) 
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At the interface between the sediment and the water 
column, i.e., at z=H(r,@), the boundary conditions are 
continuity of pressure and continuity of the particle velocity 


component normal to the interface: 


DP, U2Z-8; ea) — ee ee) (9) 


A A 
‘ (10) 
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= 
The subscripts 1 and 2 denote the water column and the 
sediment, respectively. 

The unit directional™ vector fh normal to thee peeeen 


tirethace Hs 
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where Z,f and @ are the unit directional vectors associated 


with the z, r and @ directions, respectively. In the case of 


1 


a small bottom slope, the boundary condition of Eq. (10) can 


be approximated by 
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The small slope approximation is accurate when 


OH(r,8) tee ( r, 8) 
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Or 
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Bem@lowing Eq. (9) and Eq. (12), we obtain the following 
boundary conditions at the water-sediment interface for the 


normal modes: 


1 02,,(2=H;r,8) LOZ. Zane, UO} 
fee ee eee SS ee ee (14) 
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These boundary conditions hold for each individual normal mode 
because they are orthogonal functions. 


The boundary condition for p at 270 is 


Del 2270, Ue (16) 


This implies that the ‘normal modes and their horizontal 


derivatives are also zero aS Zz > oOo: 
Z,,(Zz7~;r,9) = 0 (79 
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i O25, ae oe 19 
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C. MODE COUPLING COEFFICIENTS DERIVATION 

Substituting Eq. (2) into Eq. (1), multiplyingaia, 
Zrn(z;r,8@)/p, integrating over the entire depth and finally 
rearranging terms, we obtain the coupled mode equations 


governing the mode amplitude functions: 


[V,2+k2,,(r,0)] R,(r,8) = -}> [A,,(z, 8) V,R, (rz, 8) +B,,,(r,8) R(x, B)] 
Nn 


(20) 


where the two mode coupling coefficients are defined as 


Ae = 2f—z, (27, ON Ve an ze Onell (21) 
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and 


Bunir,98) = 


I ; 
2m (Zi 2,0) Vp, (zi 2,8) dz (22) 


V, is the horizontal gradient operator, i.e., 
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A A (23) 
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The mode coupling coefficients are measures of the 
Significance of exchange of acoustic energy between modes 
resulting from horizontal variations in the medium. As the 
variations become stronger, the coupling coefficients become 
larger and so is the energy exchange. In the case of a 
completely range independent medium, the coupling 
coefficients are identically zero and the RHS of Eq. (20) 
vanishes. In such case, the modes propagate independently of 
each other.For range-dependent cases, the neglect of mode 
coupling leads to the adiabatic approximation [Ref. 5]. 

Cylindrical spreading can be removed from the coupled mode 
equation (Eq. (20)) by replacing the mode amplitude function 


Pueem) with P (r,@)/r'*. The result"is 
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III. THE NUMERICAL MODEL AND EXAMPLE RUNS 


In this chapter, the procedure to upgrade the Chiu-Ehret 
model [Ref. 6] is discussed. The upgrade has entailed the 
derivation of alternative expressions for the mode coupling 
coefficients and the generation of new code to compute these 
coefficients based on the alternative expressions. 

The numerical results from two simple example model runs 
associated with two different bottom slopes are also presented 
in this chapter. Both cases deal with upslope propagation in 
1sospeed layers. These runs have allowed for an examination of 
the coupling between modes caused by bathymetry change. In 
addition, they have allowed for an examination of the validity 


of the adiabatic approximation. 


A. THE CHIU-EHRET APPROACH 
In the far field, i.e., kr>>1 , the coupled mode equation 


(Eq. (24)) for the mode amplitude functions, can be recast as 


0" nae 
Le a eat ke = 
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(25) 
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In the Chiu-Ehret model [Ref. 6], P, is decomposed as 


i 


C10)4- Cl aeOyme. 


ie (26) 
One, 0) = [k,(r,8) dr 
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where U, is the slowly varying complex envelope of PP, 
modulating the rapidly varying two-dimensional (2D) adiabatic 
solution, i.e., exp(j¢,), and ¢@, 1s the adiabatic phase. The 


Chiu-Ehret model numerically computes the envelopes U, instead 


of P, using Runge-Kutta schemes. 


B. ALTERNATIVE EXPRESSIONS FOR THE MODE COUPLING COEFFICIENTS 

For simpler numerical implementation, the expressions of 
the mode coupling coefficients in Eq. (21) and Eg. (22) are 
rewritten in alternative forms. These alternative forms do not 
require integrations of expressions involving the horizontal 
derivatives of normal modes. In the following, the derivation 


of these alternative forms is presented. 


1. VECTOR MODE COUPLING COEFFICIENT, A,, 


a. Case of mé£n 
Applying the horizontal gradient operator V, to 


both sides of the depth equation (Eq. (3)), we get 


ICs 


#V,Z,, 
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Multiplying Eq. (27) by Z,(z;r,8@)/p and then integrating over 


the entire depth, we get 
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In order to recast the first term of Eq. (28) ina form useful 
for this derivation, we first use integration by parts twice 
with respect to z. The resulting expression, after some 


lengthy manipulations, is 
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(29) 
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Again, the subscripts 1 and 2 denote the water column and the 


sediment, respectively. 


Application of the boundary conditions Eqs. (6), (7), (8), 
ener  (i5), (17), (18) and (19) to Eq. (29), yields 
subsequently 
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mMemlacing the first term of Eq. (28) by Eq. (30) and then 
making use of the depth equation (Eq. (3)), we obtain the 
following alternative expression for the vector mode coupling 


coefficient, for mén 


Als, 
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or equivalently, in light of the boundary conditions Eq. (14) 
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(32) 


The above expression only involves Z, and not their horizontal 
derivatives in the integrands. Therefore, the corresponding 
numerical evaluations are more efficient. 

The last two terms of Eq. (32) express the direct 


contribution of bathymetry change and sediment properties in 


= 


Ann: They were excluded in the previous model but are included 


in the latege vercrone 


NG 


b. Case of m=n 


In order to derive an expression for the vector 


mode coupling coefficient for m=n, we differentiate the 


orthonormal condition Eq. (4) using Leibniz rule. The result 


is 
oe A 
ea OD) eee nee = (33) 
1 AL 
- V,H(r,0) (—-—) 2Z,,’ 
: pee | Z= er, 0) 
MNewe that, this coupling coefficient is zero for a flat 


horizontal bottom. The latest version of the model has 


miemnuaea Chis new term. 


2. SCALAR MODE COUPLING COEFFICIENT, B.,,, 


Taking the horizontal gradient of both sides of Eq. 


oe), te. ; definition of the vector mode coupling 


coefficient, and applying the Leibniz rule ene 


differentiation of a definite integral, we get after some 


manipulations, the following expression for the scalar mode 


@oupling coefficient: 


Ly 


B_,(r,8) = = Vi Ann (£8) : fs V,Zioeaa = 
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(34) 
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There is a unique property associated with a complete set of 
orthonormal functions, called the "closure relationship." For 


the normal modes, this relationship can be expressed as 


3 ini Zi 18) ZAz';r,8) = 8(z-z") (35) 


Taking the horizontal gradient of both sides of the closure 
relationship, multiplying by Z,(z;r,@) and then integracume 
over the entire depth, we get, after some rearranging of 


Cerms, 


V,2pe? 2 0) == abe cea 2) ee) (36) 
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eepstituting now Eq. (36) in Eq. (34), we finally obtain the 


following alternative expression for B,,: 
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Eq. (38) is valid for both the m#n and m=n cases. The last 
term of Eq. (38), is new in the model. The magnitude of this 
mode coupling coefficient is generally much less than the 


magnitude of the vector coefficient. 


C. NUMERICAL IMPLEMENTATION 

The major part of the model upgrade was the replacement of 
the old routines with new ones for the computations of the 
rederived mode coupling coefficients according to Eqs. (32), 
fe eand (33). 


These new routines are contained in a program called 


"“sedbot" and are listed in Appendix C. Normal modes and the 


IBS, 


horizontal gradient vector of the wavenumber as function of 
position are required as input to "sedbot." 

The normal mode field is created by a program called 
"modes", whereas the horizontal gradients of the wavenumber 
are calculated in the program "kder." These two programs are 


listed in Appendices A and B, respectively. 


D. EXAMPLE RUNS 

For both example runs, the medium is taken to have two 
isospeed, isodensity layers separated by a constant-slope 
interface. Sound speed is taken to be 1500 m/sec in the water 
column and 3000 m/sec in the sediment. Density is taken to be 
1000 Kg/m’ in the water column and 2500 Kg/m*’ in the sediment. 
The source is taken to be harmonic in time, with a frequency 
of 100 Hz, and is located at a depth of 50 m. The coupled mode 
solution along two radii, 90° and 45°, are displayed and 


discussed. 


1. BOTTOM SLOPE = .001 RADIANS 
The bottom slope in this first case is taken to be 
.001 radians. The water depth is 100 m at the source location 
and 70 m after 30 km in the y direction (see Fig. 2). 
At the source location there are twelve trapped modes in 
the water layer. Only eight trapped modes exist at the 


location 30 km upslope. 
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Figure 2. Geometry of the first example case with a 


constant slope of .001 radians along y-axis (a), and a plane 
view showing the @ = 90° and @ = 45° propagation paths (b) 


Zi 





fx) 
Q 5 g 
= punt a= = 
= re oe 
= FP 
O 
3 4 
fr) 
A. 
© 
: 3 
rn wan avd 
©, Do, aw ~ INV ahr se vier adem eye me a 
= 
: 
= 
(a 1 
© 
= 

O 

INTEGRATION DISTANCE (KM) 
Figure 3. Envelope amplitudes of the first eight trapped 


modes in the 3D coupled mode solution along the path @ = 90° 
for a bottom slope of .001 radians 
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Figure 4. Envelope phases of the first eight trapped modes 


in the 3D coupled mode solution along the path @ = 90° fora 
bottom slope of .001 radians 
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Figure 5. Envelope amplitudes of the first eight trapped 


modes in the 3D coupled mode solution along the path @ = 45° 
for a bottom slope of .001 radians 
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Figure 6. Envelope phases of the first eight trapped modes 


in the 3D coupled mode solution along the path @ = 45° fora 
bottom slope of .001 radians 
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Fig. 3 and Fig. 4 show the amplitude and phase of the 
modulation envelope, U,, for the first eight modes travelling 
in the upslope y direction, i.e., along the path @ = 90° (see 
Fig. 2). An upslope enhancement is noticed, especially for the 
higher order modes, as they propagate into shallower water. 
The phase of the envelope, which is the phase deviation from 
the 2D adiabatic approximation, is very small (about 11° 
maximum). The amplitude fluctuations are between 15% and 30% 
for all the modes. In light of the small amplitude and phase 
Fluctuations, the adiabatic approximation can be considered 
reasonable along this propagation path. 

Fig. 5 and Fig. 6 show the amplitude and phase of the 


modulation envelope, U,, for the first eight modes, along the 


propagation path @ = 45° (see Fig. 2). Here, the upslope 
enhancement is significantly less and the fluctuations of the 
amplitude of the higher order modes at greater range are 
slightly larger than along the previous path. We speculate 
that this slight increase in the fluctuations is due to that 
more interacting modes remain trapped in the water column at 
longer ranges along this path. The higher order modes have 
large phase deviations from the 2D adiabatic phases. These 
large phase changes correspond to significant horizontal 


refraction of the wave fronts due to the existence of a 


transverse gradient in the bottom bathymetry. 
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Figure 7. Geometry of the second example case with a 
constant slope of .002 radians along y-axis (a), and a plane 
view showing the @ = 90° and @ = 45° propagation paths (b) 
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Figure 8. Envelope amplitudes of the first eight trapped 


modes in the 3D coupled mode solution along the path G-= 90° 
for a bottom slope of .002 radians 
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Figure 9. Envelope phases of the first eight trapped modes 
in the 3D coupled mode solution along the path @ = 90° fora 
bottom slope of .002 radians 
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Figure 10. Envelope amplitudes of the first eight trapped 


modes in the 3D coupled mode solution along the path @ = 45° 
for a bottom slope of .002 radians 
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Figure 11. Envelope phases of the first eight trapped 


modes in the 3D coupled mode solution along the path @ = 45° 
for a bottom slope of .002 radians 
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2. BOTTOM SLOPE = .002 RADIANS 

For this case, the same isospeed, isodensity, wedge 
shape waveguide is used, except the bottom slope is now 
doubled (.002 radians). The bottom depth at the source 
location is now 100 m and shoals to 40 m after 30 km away from 
the source in the y direction (see Fig. 7). At the source 
position there are twelve trapped modes, but in 30 km upslope, 
there are only five trapped modes in the water column. 

Fig. 8 and Fig. 9 show the amplitudes and phases of the 
modulation envelope, U,, for the first eight trapped modes 
travelling upslope in the y-axis direction, i.e. along the 
path #9 = 90° (see Fig. 9). Upslope enhancement is much 
stronger than the previous case, especially for the higher 
modes. The fluctuations in amplitude is about 50% in some 
modes and in phase more than 20°. Thus, the adiabatic 
approximation would induce considerably larger errors than the 
case of a .001 bottom slope. 

Fig. 10 and Fig. 11 show the amplitudes and phases of the 
modulation envelope, U,, for the same eight modes along the 
propagation path, i.e., @ = 45° (see Fig. 7). The horizontal 
refraction phenomenon is much stronger here than for the case 
of a .001 slope. Along this path, the adiabatic approximation 
would also induce large errors. Typical percentages of 


amplitude fluctuations are about 50% for the second mode and 


a2 


30% for the third and fourth modes. The phase deviation, 


especially for the higher order modes, is also large. 
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Vis CONCLUSIONS AND RECOMMENDATIONS 

A 3D coupled normal mode model for sound propagation in 
shallow water with irregular bottom bathymetry, is developed 
in this thesis. This model can be used to examine underwater 
sound propagation involving significant bottom interact 107maeen 
this model, sound speed is allow to vary in three dimensions 
and water depth and sediment properties in horizontal 
Lecat Hem. 

It 1s shown here that, for a frequency of 100 Hz, the 
adiabatic approximation is valid only for very mild bottom 
slopes. Typical errors for a slope of .001 radians are 15% in 
mode amplitude and 10° in its phase. For a slope of .002 
radians, the errors are significantly larger. 

The model presented in this thesis 1s capable of 
Simulating the interactions of the normal modes as they 
propagate in complex environments. Propagation phenomena like 
mode-mode interaction, horizontal refraction and slope 
enhancement can be examined using this improved model. 

In the development of the present model an approximation 
(Eq. 12) in the bottom boundary conditions is used.) tne 
validity of this approximation requires that the slope must be 
much smaller than unity (Eq. 13). In order to be able vee 
handle very steep bottom slopes, 1.e., order one slopes, one 


needs to use the exact form of the bottom boundary conditions 


34 


(Eq. 10). This could make the formulation more complicated but 
it should be tractable. 

Another future improvement to the model will be to include 
sound energy absorption (attenuation). One way to do this is 
by wbialic -(efelolelaliale; imaginary parts nlgig! the elgenvalues 
(wavenumbers). Lastly, a test for the accuracy of the improved 
model is needed. This can be achieved by comparing the results 


generated by this model with some exact analytic solutions. 


Sus 


APPENDIX A. FORTRAN ROUTINES FOR COMPUTING NORMAL MODES 
FIELD 
The following program creates the normal mode field, 
"amode.dat", for a given geographical area of the ocean. Sound 
speed, density, and bottom depth, are defined for every grid 
point. Given sound speed field and density, the normal modes 


are calculated using a standard mode solver routine. 


KEK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KK KKK KKK KKK KKK KKK KKK KKK KK 


* * 
x This program computes the normal modes ina given * 
ts emaear. bs 
* * 


KKK KKK KKK KKK KK KK KKK KKK KKK KK KKK KKK KKK KEKE KEK KKK KKEK KKK KKK KK KKK KKK 


c 

@ INPUT/ARGUMENTS 

Cc 1g frequency, Hz 

Cc xmax maximum position in x direction, meters 

e ymax maximum position in y direction, meters 

el n(x; ase depth at x,y position, meters 

G nx number of stations in x direction 

G ny number of stations in y direction 

e Tz, number of stations in z direction 

eS DG) initial x position, meters 

é yo initial y position, meters 

q dz step size in z direction, meters 

C Gi (nxeny ane) sound speed field in every x,y position, 
Cc m/s 
e GE (scan. tas) density field in every x,y position, 

C kgr/m*3 
C isw switch index : 1 write / 0 deo not wee 

e mm maximum number of allowable trapped 

e modes 
c 

e OUTPUT /ARGUMENTS 

c for each horizontal station x,y: 

‘@ KSqQ isi) Squared eigenvalues for each trapping 

e mode (real) 
Ee efunm Binz ,az eigenfunctions (real) 

Cc h(nx,ny) water depth (meters) 

‘o eit) sound speed profile in a specific grid 


BIG 


OL ORIG Olas Gaal Ge i@ ae 


2) 


position 


Gov) density profile in a specific grid 
position 

Ww source angular frequency (rad/sec) 

se, dy 7dz step size in x,y and z directions (meters) 

Hee Ty, 11 number of stations in x,y and z directions 
(metrers) 


program modes 


oe 30000.d0, ymax=30000.d0,nx=11,ny=11, 
nz=100,mm= 20) 

implicit real*8 (a-h,o-z) 

median ek (doenlye a2) 4c VIZ) AE (ix, my,mz) ,d (mz) 

neal*8 hinkx,ny),ksq r(nz),ksq 1(nz) 

real*8 efun r(nz,nz),efun_i(nz,nz),efun(nz,mm) 

real*8 ks(mm),x(nx),y(ny),z(nz),ksed,kwat 

logical ex 


data isw /1/ 


inquire (file='amode.dat’ ,exist=ex) 
if (ex) then 
open (unit=13, file=’amode.dat’ ,status=’old’) 
close (13,status=’delete’ ) 
endif 
open (unit=4,file=’ /home/noise/sagos/modes/amode.dat’, 
* form=’unformatted’,status=’ new’ ) 
inquire (file=’mode.sys’ ,exist=ex) 
if (ex) then 
open (unit=13, file=’mode.sys’,status=’old’) 
close(13,status=’delete’ ) 
endif 
open (unit=6, file=’ /home/noise/sagos/modes/mode.sys’ 
* form=’ formatted’ ,status=’ new’ ) 


write(6,*)’output field’ 
dx=xmax/dfloat (nx-1) 
dy=ymax/dfloat (ny-1) 

ola sal) 
pi=4.d0*datan(1.d0) 
£=224 .d0 

WZ Selo aoe 


cali wemee, tf, h mx,ny,nz,dx,dy,az,x0, vO, *,V;Z) 
nm=nz-2 
write(4) w,dx,dy,dz 


write(4) nx,ny,nz 


1f (1isw.eq.1) then 


ae) 


OO) CONG 


eG 


write (6,*)’daz= ',dzZ se meters. 

write (6,4) mx= 79) Cry — a le 2 
write(6,*)’interface depth’ 

write(6,*)’h= ',h(i,1) ,’ meters’ 


wWrite(6,*) (Chi az) iz eee 

write(6,*)’density profile, kgr/m*3’ 

write (6,*) (adi (41 22) ae 
endive 


( 
( 
( 
( 
write(6,*)’sound speed profile, m/s’ 
( 
( 
( 


Meter —o 


Ger 1 x= nen 
dew y— lL any, 
Teoumtter—@ 
ichk=0 
1f (1xte@. wand y7 cay arene 
write(4) h(1ix,1y) 
doy 2Z= 1 ne 
GGz t= Cis (eee 
Gl (352) —6 Base yee) 
enddo 
write(4) c 
write(4) d 


write (6,*)’ix=",ix,” dys y seme Genelia 
call mode(£,;nz,dz,c,djnm, ksqur,ksq 
_ Cfunw Yr, Chun) jyameliic) 


choose only the trapped modes 


cs : sound speed in the sediment (constant) 
cw : sound spedd in the water column,next to the 
interface 


CS=c (ink (h(i ay dizi) 
Cw=c (int (a (ise) diz) 


set zeros in the eigenvalues-eigenfunctions arrays 


do i=1,mm 
KS (a) =07 a0 
Glo} y= 1. al 
ens bbaiig|; st) —Oacle 
enddo 
enddo 


do i=1,nm 
ksed=(w/cs) **2 
kwat= (w/cw) **2 
if (ksed..1lt.ksq. ri) amd? heaq san) ee waits 
* then 
Lcounter=icounter+1 


sie: 


ks (icounter) =ksq_r(1) 
dow Z=1,nz 
efun(iz,icounter) =efun r(iz,i) 


enddo 
endif 
enddo 
‘@ 
if (isw.eq.1) Bee 
Wieiiee(G , ~)eel— "Suhel ) 
yaita(G «cl mire for trapped modes : ’ 
write(6,*)kwat,ksed 
write(6,*)’icounter= ’,icounter 
WindiGemige® )? Iara 2” 
Widteer( Grete (isn 1) Pik mn) 
endif 
‘eo 
Pf mee 3.end. 1y..eq-3)  Lhen 
write(6,*)’efun(iz,18) ’ 
Weite (6. *) (e@fun (12, 13) ,1z=1,nz) 
endif 
@ 
write(4) icounter 
lh | LeCOUNnEGCHr-ge.mm.and.1Counter.gt.mcnetr) 
* mentr=1icounter 
write(4) ks 
write(4) efun 
enddo 
enddo 
@ 
Pie (menerene.0) write (6.1001) mentr 
1001 format(i3,’ trapped modes, exceeds limit, increase mm 
‘e and rerun’ ) 
close (4) 
close (6) 
end 


KKK KEKE KEK KEKE KEKE KK KKK KK KKK KKK KKK KKK KKK KKK KKK KEK KKK KKK KKK KKK KKK KKK 


* * 
* The following program provides an example data input. * 
* * 


KEKE KEK KKK KKK KEKE KKK KEKE KEKE KEKE KEKE KK KKK KEKE KEKE KKK KKK KKK KKK KKK KKK KKK KK KKK 


S 

‘@ INPUT/ARGUMENTS 

e nx number of stations in x direction 
Ss ny number of stations in y direction 
e az, number of stations in z direction 
e dx step size in x direction, meters 
a! dy step size in y direction, meters 
eC az step size in z direction, meters 


ao 


CR? ae ie oes One oak ewe. 


Q 


Q 


Q 


xO initial x position, meters 


yo initial y position, meters 
OUTPUT /ARGUMENTS 
CE (nx nym) sound speed field in every x,y position, 
m/s 
Ge ( Nocany., mize density field in every x,y position, 
kgi/ ties 
(nec, ay) interface depth, meters 


subroutine data(cfi,dfi,h,nx,ny,nz,dx,dy,dz,x0, yo, «= aaa 


Iimplicie ieeal *3 saa oe 
real*8 cf (nxny, nz), dE (i) Dy nz eee 
real*8 x(nx); y (ny) zz) 


Ghe) toke= Neb 
GO many Iya, 
x (35¢) =xo+dsc*dt logins 15c— 1 
y (iy) =yotdy*dfloat (iy-1) 


bathymetry field 
h(ix, iy) =-.0005d0*x(ix) +100.d0 
sound speed and density fields 


do iz=1,nzZ 
Z (1Z) =dz*dt loa az a) 
1£ (z(iz).le.h (ase 1y eaeben 
GE (ix 1 y ,12) =. OOOS00 eGee)n- 
* .1d0*z (1z) +1490.d@ 
GE (156, 1 52) Or a0 
else 
GCE (1x ay az) = Oe a) 
Gi( 1x7 iy Zz) = 210 COM 
endif 
enddo 
enddo 
enddo 


return 
end 
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APPENDIX B. FORTRAN ROUTINE FOR COMPUTING WAVENUMBER 


DERIVATIVES 


This program inputs from "amode.dat" as created by the 


previous program "modes". Computes the horizontal derivatives 


of the total wavenumber k, at every position of the acoustic 


field. The derivative calculation requires definition of a 


computational domain. Output is to the file "kder.dat". 


KEKE KKKEK KKK KKK KKK KEK KEKE KKK KKK KEK KKK KKK KKK KKK KKK KKKEKKKKKKKKE KKK KKK 


* 
* This program assigns the source position relative to 
* the input field via ixorig and iyorig. 
* This program also specifies the radial increment for the 
* spline definition, the number of intervals, and the 
* angular increment between integration paths (dr,nr,da) 
* 
* Procedure: xy-Spline at each depth 
* evaluate dk/dx,dk/dy 
* transform into dk/dr,dk/da 
* 
KERR KKKKEKR KKK KEK KKK KKK KEKE KEKE KKK KE KKK KKK KKK KEKE KKK KKK KKK KKK KK KK KK 
€ 
pregqramn Kader 
Cc 
parameter (nx=11,ny=11,nz=100, nm=20, ndum=nzZ*nmn, 
* 1xorig=1,ilyorig=3,nwk=2*ny*nx+2*max0 (nx,ny) ) 
implicit real*8 (a-h,o-z) 
Bea os wie oe, yi 12,2) ax, Ny) ,.Koxy (6), ke(2,nx,2,ny) 
mew Ou x (ei y (iy), ang (nx, my). c(nzZ) Clink, ny,1z) 
real*8 wk(nwk) ,efun(ndum) , von (cae ) 
character*20 filename 
logical ex 
Gc 
‘Ge —--*- ef --e- eee ewe ewnewe ewe ewe he we ew ew ee ee we ee ewe ee ee ee ew ew ew ee 
© Open statements 
e a 
@ 


inquire (file=’kder.dat’,exist=ex) 
if (ex) then 


open (unit=13, file='kder.dat’,status=’old’ 


close (13,status=’delete’ ) 
endif 
inquire (file=’kder.sys’ ,exist=ex) 
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) 


* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 
* 


Lf (Cex)) cE nen 
open (unit=13, FILE=’kder.sys’,status=’old’ ) 
close(13,status=’delete’ ) 

endif 

inquire (file=’efun orig.dat’,exist=ex) 

if (ex) then 
open (unit=13,file=’efun orig.dat’,status='old’) 
close(13,status=’delete’ ) 


endif 
open (unit=4,file=’amode.dat’,form=’unformatted’, 
* status=’ old’ ,err=2 00m) 
open (unit=14,file=’kder.dat’,form=’unformatted’, 
* status='’new’ ,err=2002) 
open (unit=24, file=’efun_orig.dat’ , form= 
* ‘unformatted’ ,status=’new’ ,err=2003) 
open (unit=6, file=’ kder.sys‘* , £oum— Ecrmaceeas 
* status=’new’ ,err=2004) 
@ 
c input w (rad/sec), dx,dy (meters), dz (meters) 
read(4) w,dx,dy,dz 
Cc 
c input: number of x indices, no. y indices 


e number of modes, TOTAL vertical increments 
read (4) Gixxe emy ey Zit 


é 
write(6,1009) w 
write(6,*)’ dx (eddy (m)7ezin) 
WELte (6,* )@ doc, dy ,dz 
1—£ (Mine nm) welte (6, *) m= 2m, em— am 
if (nxX.nesnxt) write (654) tx 2 oe, ees ee 
if (ny.ne.nyt) write(6,*) ny=' ,ny,°  nyt= ave 
if (nzpl.ne.nzt) write(6, *) nzpol=" nzpl,’ _nzt= sae 
a 
pi = dacos =iad0) 
@ 
C =--------- distribution parameters for spline------------ 
dda =a cad 
da” = dda sey 1402 


di = 20007a0 
c number of points in spline 

ew (re-do ) tebe cia | 

write(6,*) nr,’ Spline locations with interval— [dm 
c number of radial paths 


tang = datan( .5*dy*(ny-1) -/ 4dx* G1, 
e na = 2*idint(tang/da) + 1 

na=3 

write(6,*) na,’ integration paths for da=’,dda,’ deg’ 
Ge ee oa eae ee ae a ae ew we ee So Se eS ee Se eS ee ee eee ee ee ee ee ees aes ce ee es cee ce ee cee ee ce ee ef ee 
@ 


c horizontal field grid in meters 
dior 11° as<= 19 aes 
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Q 


Q (2) 12) elie 


are 


ASA aso 


Ji x(ix) = dfloat (1x-ixorig) *dx 
ieee dy — ey any: 

eevee — at Lode (ly 1yorig ) *dy 
VWHeniomeGn. Se xa eicie: (yd) 4 x (TX) 
WieinentomemiecesmacdmGe =: (" yids) § -' eyviny), 


read in c-field and 
calculate angle (ccw from x-axis) 
OMG SE Sloe AR a>-< 
do 14 iy=1,ny 
full sound speed profile(nzpi1) (0-5000m) 


read(4) h 
read(4) c 
read(4) dens 
read(4) icounter 
read(4) hork 
read(4) efun 
create file to obtain initial conditions 


if (ly.eq.iyorig.and.ix.eq.ixorig) write(24)efun 
do 1Z=1,nzZ 

eee 17) = Co (1) 

enddo 

hee x 0 1 xOri.g) then 


eevee Ge vi@mie |. ONG (ix ay) = 1/2. 
teianweliaaivorig) rang (ix, iy). =piy 2. 
else 

ang(ix,iy) = datan((y(iy))/(x(ix))) 
endif 


14 continue 


calculate derivatives from first station below surface to 
DOtELom 


do, 1007 7Z=1 nz 
do 110 ix=1,nx 
Gordie yi ny 
wavenumber is in rad/m 
Lake imigiec inven = W/OL (1x, 1V, 27) 
fit bi-cubic spline to iz-th level waveno. 
ic = nx 
Coll TOCeGeU(K cnx ,y Ny ce,C, why 1er) 
1f (ier.ne.0) write(6,1001) ier 
use spline to evaluate cartesian derivatives at each grid 


point 
transform derivatives into cylindrical coordinates 
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tad 
20 


A916) 


LOO: 
1OOZ 
1003 
1004 
LOO 5 
ALONONS 
gO le 
LOO 
ION One 


* 


do 1205 1—1) nx 
Gow way 1 ty, 


call dbcevl (x, nx jammy (Ke, te) Gee aa, 


a 


Call eyitiang (ce 


,Kdaxy jee 
(ler.ne.0) 

write(6,1007) 151), x(a 1er 
Aiobevell 2) 


Kea cheaey 74.2 


OL 
PSD y1ceh)) 
Fe OU ale cree al van ZN) 


continue 
continue 


ae 


(1z.eq.1) 
write(6,*)’ 


then 
For first subsurface layer’ 


write(6,1005 


write(6,1006 


((ct (1x,1y,1) ,1x=]=1,nx) joy ee 


write(6,1011 


) 
) 
) 
write (6,1016) 
write(6,1003) 
write (6,1002) 
write ) 
write ) 
write 


endif 


continue 


write(14 
write(14 
write(6, 
write(14 
write 
write 
write 
write 
write 


(6 
(6, 
( 
( 
( 
write ( 
( 
( 
( 
( 


6, 
Sy 
OF 
Sy 
oF 
6, 


write 
write 
write(6é, 
write(6, 
write(6, 


Goto 202 


format ( 
format ( 
format ( 
format ( 
format ( 
format ( 
format ( 
format ( 
format ( 


((k(1x,1y) (1 x= $ise) eye 


((kd(ix,1y,1,1) , lxX=1) mx) fay ay 
(6,1004 
(6, 1002 
( 


Sey 


((Kdigoex ay 1, 2) 


pix=1 nk) | 1y=ny , 


DONS aS FL ZOLG ALG 
da,na,dr,nr 

"da ma, dr ; nn 
kd 


) 
) 
* ) ada naar ,nr 
) 

*)’ For bottom level’ 
1005) 
1006) 
LOnLd } 
1016) 
LOO 2 
Loo2) 
1004) 
1002) 
* } 


1008) 


((cE (1x, 1y, nz) Fie), mse y =e 


({(k(1x, ly) 7 ixa1 na — ee) 


((kd (ix,ay,nz, 1 as se) ty a 


((kd (1x,1Y,nZ, 2) 74 xk=1) oe yan ee ee 


(kd(4,3,12Z,1) .22=1,n2, 2) eas 3 nz 


Ler: * 38, “for abcecu, xsy-spline®) 
(ix 61275 )9 
Gic/ Ginna 

ic) Glau 
c (m/s)’‘) 
(2x,£8.3)) 
(2x,£8.5)) 
ix, ly,X, Y¥> ler fom dbcevwr: 
at 1x,1y=4,3 Gk adnizy 7 linus 


i anes 
jee) 


ee 
(rad/m) 


OTST 2h Ge ie 3) 
(1x dli.s5 7) 
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moo > format (” £r er ae rad/sec :',d14.7) 
mono tommat(’ k’/7 (4(3x. — 4) i 3 ape Gd 74) /) 
mOdd format (’ k Aone : 
C 
ee ee Pee oe Ee le 
c close statements 
eee SESS 2 6 eee nee nee ences 
C 
2001 filename=’amode.dat’ 
GOeoeZ0 10 
2002 filename=’kder.dat’ 
Gocco Zou 
2003 filename=’efun orig.dat’ 
GOE@O °2010 
2004 filename=’kder.sys’ 
2010 write(*,2011) filename 
2011 format(’ ERROR OPENING FILE ’ A) 
2020 close(4) 
close (14) 
close (6) 
end 


Subroutine cyl(ang,x,y,r,a) 


c polar transformation subroutine 


implicit real*8 (a-h,o-zZ) 

r = x*dcos(ang) + y*dsin(ang) 
a = x*dsin(ang) - y*dcos (ang) 
return 

end 


45 


) 


APPENDIX C. FORTRAN ROUTINES FOR COMPUTING THE COUPLING 
COEFFICIENTS 


RRR KEKE KEKE KEKE KKK KEKE KKK KKK KEKE KEKE KEKE KK KKK KK KKK 


This program manages the subroutines "Subl1.f", "sub2.f", 
"sub3.£", "partial.f" which compute the two mode 
coupling coefficients. The input is from "amode.dat" 

and "kder.dat", specifically the modes, horizontal, 
wavenumber, horizontal derivatives of total wavenumber, 
bathymetry and density. The output file is "mcoupl.dat". 


> OO 


KEKE KEK KKK KKK KKK KKK KEKE KEKE KKK KEKE KEKE KKK KKK KKK KKK KK KKK KKK KK KEKE 


rhol : water column density-constant in depth (kg/m‘3) 
rho2 : sediment density-constant in depth (kg/m‘*3) 
dz : vertical step size (m) 
nm : maximum number of trapped modes,in the water column 
nx,ny : number of stations in x and y directions 
h : bottom bathymetry 
zb : acoustic pressure eigenfunctions,at the interface 
depth 
zobm1l : acoustic pressure eigenfuctions,one step size 
above the interface depth 
zbm : zb for the mth mode 
zomm1l : zbml1l for the mth mode 
Zon > zb Eor the neh mode 
Zpnmi >) Zomietor Eheait a emede 
ar : range component of the first coupling coeff. 
aa : angle component of the first coupling coeff. 
cr : correction at the range component of the first 
coupling coemms 
ca : correction at the angle component of the first 
coupling coeff. 
b : second coupling coeff. 
k +: square of horizontal wavenumbers (eigenvalues) of 
the modes 


GG ey i Or Or Ore ae One ey (Oo Or a OG OO). "Or OP ae be ee oh ee oe oe 


program sedbot 


Q 


parameter (nx=5,ny=5,nz=100, nw=2*nx*ny+2*max(nx,ny) , 
* nm=20) 
implicit real*8 (a-h,o-z) 
real*8 hinx,ny), zom( nx; ny) bin mm moe ay) em merge 
real*8 zbmml (mx, ny) , Cr Gall am, nx, ny) caunm, an, wee 
real*8 ar(nm,nm,nx,ny) ,aa(nm,nm,nx,ny) ,Kn(nx,ny) 
real*8 zbn(nx, ny) 7, zZbnmipG@ix me (2 ee 
real*8 cri (nx, ny) , call (1c) my ys se)  yay) 
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Q 


Beal seme (ioe hy) aio iloa | Wea yeeeZiont Tim, Moc; ny) 

Beal omelet oc pve | 2onovn nec my je ior (mx eny) , hoal(nx, ny) 
ieeane* So millon (oc) | inioyn oe, ye odeiiasc, iy), Ki (11m) 

Bedi wom molto, NY) ,bNO2 (Nic my ki nx, ny, Zz) 

hei omesn ize Gd (az), erin (naz, nm), kedimx, nys;nz. 2) 

integer icounter (nx,ny) 

logical ex 


inquire (file=’mcoupl.dat’ ,exist=ex) 

if (ex) then 
Openiare=13 file=’meoupl dat’ ,status='old’ ) 
close (13,status=’ delete’ ) 

endif 

open (unit=4, file=’ /home/noise/sagos/modes/amode.dat’, 

* form=’ unformatted’ ,status=’old’ ) 

inquire (file=’ coupl.sys’ ,exist=ex) 

if (ex) then 
Geen (Unit=13 > file="coupl2syvs*,status='old’ ) 
close (13,status=’delete’ ) 


endif 

open (unit=6,file=’ /home/noise/sagos/modes/coupl.sys’, 

* form=’ formatted’ ,status=’ new’ ) 

open (unit=8, file=’ /home/noise/sagos/modes/mcoupl.dat’, 
* form=’unformatted’ ,status=’ new’ ) 

open (unit=14, file=’ /home/noise/sagos/modes/kder.dat’, 

* form=’unformatted’ ,status=’old’ ) 


read(4) w,dx,dy,dz 
rewind 4 


Ay ExXOrlG VyOri1Lg 
4) da,na,dr,nr 

write(8) ixorig,iyorig 
gyadar Na,adr,nr 
4) kd 


do. jt <= 1. 
K(1x) deer df loat (1x-1x0rig) 


enddo 
Come y— lay 
y (iy) =dy*dfloat (iy-iyorig) 
enddo 
Witmesiom ns x ramge. (xd) yp, § (ax) ." ) 
Wihitentone) “vecange: “(7 jyil)! .)oyviny) ,!)' 


computation of first mode coupling coefficients 


4‘7 


eee 7 ee 


do sam 

den m= pain 
read(4) w,dx,dy,dz 
read(4) nxx,nyy,nzz 


elo eal ee. 
don ivy 
read(4) h(ix,1y) 


1, is the last station in the water column 


T=int (hh (isc yy ez 
read(4) cs 


calculate the total wavenumber 


Glia) aba il eee 
kK (1x ly paw eee 
enddo 
read(4) d 
va) gv@yal (alae abot a.) 
ChoOZ (ioe y ) =a) 
read(4) icounter(ix,iy) 
read(4) kk 
Kem doce) clear) 
| aml Galp.eeeal. 1c eal) 
read(4) efun 


if (m.eq.n) then 
aie (il -sineedecennye) = Oral) 
aal (Mii de days) Oral) 
else 


trapezoid integration, to find the integral part of 
first coefficient 


sumx=0.d0 
sumy=0.d0 
denom=km (5c) 57a ken (ise 


do iz=1,nz 
a=K{ix,1V,12) *ebuni iz a) 
* efuniaz, im) aie) 
sumx=Sumx+a*kd(ix,iy,1iz,1) 
sumy=Sumy+a*kd(ix,ly,12Z,2) 
enddo 


1f (m.eq. Tl72and nm veg aiG) 

* WELte (6, *) > tea Psunis, Sun ya ec lar 

“ Sumx, Sumy 
ar(m,n,ix,iy) =4.d0*sumx*dz/denom 
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OO 1 


aa(m,n,ix,1iy) =4.d0*sumy*dz/denom 
iMedia) - and. .eq lo ) 


x Wiese Gees Siac ey woe (hy eelO, 1x, ay)! , 
* ieee (ls Lo, 1x, ry) 


endif 


Along || aloe alla) =i 6 ol Gale cal) 
ZO ea) = LU (a Lm) 
Zio aed) =e eta (3 Tay) 
Zone (ax iy) =erun{(a-1. nn) 
enddo 
enddo 


Sawaal, nec mye x,y pipe, pas, mv, npc, apy Cc, wit) 


Gai paredal (Zon, 1x, ny, <<, y, Zonor, Zolnpa, lw, 


* ZoneX. Zonpy , Cc) wx) 


if (m.eq.17.and.n.eq.16) then 
WElte to. *) abn = ’ 
mocimine ( 6 MOO 1 Tez) ,1x=-1,nx) , 


* I= ee) 


Miele Gee) Zao a=" 
Vigiee| Gn HOG) (clmpr (ax, LY) > i x=1 50x), 


© al 4901) pe Oa ly) 


write(6,*)’zbnpa =’ 
ViebEe Oo, OO) ( (Zonpna (7,17) ,ix=L, nx), 


* uy 10) nelly lle 


POmMiee( 52x, el2.4) ) 
endif 


tiem iemimimmeat ll Sull(rhol, vidoe ,dz,zom, zon, zommi, 
* Zone reie ea kit, ciel Mm eZoOmMon Zon pa ,.<, 7) 


iPod meth he sitoZ (anol nr MO2 «2.0m. Cr l,-Cal, 
x Tein ot pal, 2c, Y ) 
Gigmaia<= Ie mac 
domiy= 1 ny 
he ane ey =O rds 21) 
Cavan dX, iy) Ca L (scp a) 
enddo 
enddo 


Pin Med hy and. te. lo) then 
write(6,*)’checking quadrature :’ 
Wise (Git) elie (Joey) =” 
Vio ommlOOmeca Misty eG doo, iy) yd x=), 5) i y=1, 5) 
Viator si aati, 16,1, iy) =" 
Viti ommEOM \iaee (1) pho , ix, ) , 1 x=1,5), 1y=1, 5) 
endif 
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100 


Gone se sees 
GO se — 1 ny, 
ar(M,n,1ix,1y)=ar(m,n, ix, ly) er m,n eee 
aa{(m,n,1ix,1y)=aa(m,n, 1xX,1y) 7ea (m,n, 2 ee 
enddo 
enddo 


rewind 4 
enddo 


Gi@e disc= eine 
Glo aby 11, sal 
yale (idly slove alhi2) —valeyol | Lee aalNg 
enddo 
enddo 


enddo 
computation of second mode coupling coefficients 


do n=1,nm 
do m=1,nm 


call sub3 (rhol,rno2, 2b) abl ar, aa,m, ml, ears 
* zbnpr, zbnpa,x, y) 


do ix=1,nx 
dey =aea 
If AaCOUNE CR 156 ia) melee yee. 
* 1counter(ix,iy).1lt.m) emem 
lech aly pipe 2h) 0) LG 
else 
ls Givala alec, aly) Slo abe aie, 
endif 
enddo 
enddo 


enddo 
enddo 


write(6,*)’checking the uejels coupling coefficients :’ 
write (6,*) an(16)17 72, — 

write (6, 100) (lant ee ae iy), ixen,5) ,dy=2,5) 
write(6,*)’aa(16,17,ix, iy) = 

write (6,100)%(aa(16,17 fase gO) ek Ab paar ah, 3 

write (6, *)” ble gas iy) = 

write (6,100) ((b(16, 17 jae ey ee 
format (5 (1lx7etZes 
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write(8) ar 
write(8) aa 
write(8) b 


close (4) 
close (8) 
close (14) 
end 


oval 
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Q + + + + * % 


* 
This subroutine computes the vector mode coupling * 
coefficient correction, due to small bathymetry * 
variations between two different modes (m different * 
than n. rt 

* 

* 


KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKKK KKK KEkK 


subroutine subi (rhol, rho2 ,dz, zbm, Zon, zon Zeume 
* erl,cal,km,kn,m,n, zbnpr pzbmparee a 


parameter (nx=11,ny=11,nw=2*nx*ny+2*max (nx,ny) ) 
implicit real*8 (a-h,o-z) 

real*8 ¢(2,nx,2,ny) ,wk (mw) cri (moe ny) cal mec a 

* zbm(nx,ny) 
real*8 zbmmi (nx,ny) , zbn (nx,ny) , zbnmil (nx, ny) , kn (napa 
real*8 zbnpz (nx, ny) , ZompzZ (nx) oy) | 2Zbeeap Se le nee 


* km (nx, ny) 
real*8 zbnpzpa (nx, ny) , zbnpzpx (nx, ny) , zonpzey (nx, 
x zobnpr (nx eis 


real*8 zbnpa (nx, ny), chol (ax) ay)yenoZ nsw 
real*8 x(nx),y (ny) 


do 1x=1,nx 
Cp sales 1k ayy 
ZbnpzZ (1x, ly) =(zbn (x, ty) Zone (scan) az 
ZompzZ (ix iy) = (Zom Gis, ty) Zon Gisele ol 


enddo 
enddo 
call partial (zbnpz nx, ny 267 )2bapzp, zoupzea ms 
* Zobnpzpx, ZbnpZpy , Cys 


if (m.eq.2.and.n.eq.18) then 
WELECG( 6G) *) Zone z (3) 3) —. , Zbnpz (3 , 3) 
WEIte (67%) zbnpzpre( 3a) — 7 aomwzor (se) 
write(6, ee. 3)=’, zbnpzpa (3 3) 
(6 f 
(6 , 


write *) Zope (47s) eet, 3) 
write *)* zonpal3 7 3) — pZioni@anesers)) 
endif 


do ix=1,nx 
Ge iy = ileal, 
F=dSqre (x (90) 42 yay) <2) 
tf (role ed 20 cere m00 
Crl (1x,4y) =ZbnpzZpin (sc, ya ee Zions oc minyees 


* (1 .d0-rho2 (1x, 47) /mnol4i ix, 1 yoy enol (ix) tye 
* ZoMmpZ (ix, dy) *zbnpr (xy 
* (ih G.Oy/aehi@de (aloes ee =), d@/ehno2 (isan) 
Cal (ey) — ion eee ly) *zbnpzpa (ix,1y) * 
* (I aO= nhO? (156, v7) 7 sinome ele ey aie, 


eZ 


* (72N@ Wie) 1s eZ In Zieis, Ly) * 
ZONWa Uy joe 
x Pie Gey aiolmiecweky toe Gley tenO2 (13, dey) .)./ 1 
Gib cy erste toe) ) ely ecm, ly) —km(ix, iy) ) 
Gab ioe, Ay) ea lea day) eco (isa, 1y) -km (ix, iy) ) 
100 continue 
enddo 
enddo 


* 


return 
end 


oe. 
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* * 
* This subroutine computes the vector mode coupling * 
* coefficient correction due to small bathymetry * 
* variations, in the case of m equals n. * 
* * 
KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KKK KEKE KEK KKK KKK KKK KKK KKK KR KKK 
e 
subroutine sub2 (hol, eho2, Zbm, a, Gi ieeade Tiberias 
* hpr,hpa, xne 
‘@ 
implicit real*8 (a-h,o,zZ) 
real*8 cri (nx,ny) ,cal (nx;ny) ,zbmi(nx ny) , rho (isc oe 
* rho2 (nxF7ies 
real*8 hpr(nx,ny),hpa(nx,ny) ,x(nx),y(ny) 
Cc 
Com aie 
Glee ayy ib ich! 
r=dsqrt (x(1x)**2 + y(iy) **2) 
if “(erie d=20) egerenlMHE 
Cri(ix, iy j= (ede) wae ee ay, 
* -1.d0/rho2 (i1x,iy) ) *(zom(ix, iy) **2) *fapr Cee 
cal (ix,1ly)=-(1.d0/rhoil (ix,1iy) -1.d0/rho2 (1x) 250m 
* (Zion (ise dhy |) 72) ioe isc nya ae 
100 continue 
enddo 
enddo 
é 
return 
end 
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* 


* This subroutine computes the scalar mode coupling 
* coefficient small bathymetry changes included. 


* 


* 
* 
* 
* 


KKK KKK KKK KKK KKK KKK KK KKK KKK KKK KK KKK KEKE KKK KKK KKK KEK KKK KKK KKK KEKE 


é 
subroutine sub3(rhol,rho2,zb,h,bl,ar,aa,n,m,hpr,hpa, 
* Zim, Zonpa, , y) 
e 
parameter (nx=11,ny=11,nw=2*nx*ny+2*max(nx,ny) ,nm=20) 
iin etedieeneal* 3” (ean, OZ) 
mocks Nal? ee, 2, lly), Wik (NW); Zo (am, nx, mya, hinx, ny), 

a bits iy) 
real*8 hpr(nx,ny),hpa(nx,ny),ar(nm,nm,nx,ny), 

* aa(nm,nm,nx,ny) 
real*8 zbonpr(nx,ny), zbnpa(nx,ny) ,ern(nm,nx,ny) 
real*8 ean(nm,nx,ny),sum(nx,ny),armn(nx,ny), 

* aamn (nx,ny) 
real*8 arpr(nx,ny),arpa(nx,ny) ,arpx(nx,ny) ,arpy (nx,ny) 
real*8 aapr(nx,ny) ,aapa(nx,ny) ,aapx(nx,ny) ,aapy (nx,ny) 
real*8 x(nx),y(ny),erm(nm,nx, ny) ,eam(nm, nx, ny) 
real*8 Eee ae MiloZ ee, day } 

! 
mo 7 
mp=16 
G 
if (ed eps ote ae Sy then 
write(6,*i@ rhel* 
ReUaG, pr a a Et Yale ee LYS 
WELbeE (6, *)° hel” 
aes (G. Eee NO ee et) 
Wirieets,*) npr’ 
WEEN. oa asia ly), 5 Bp, IRIN rma ba’ ts eg ed) 
Wiebe GO. Io Zoe" : 
eS eQUIS, 100) (8 b (my, ise a es 
write (6 were are 
wr itete. ee loyal i, iy), ix=1,5),iy=5,1,-1) 
endif 
C 
do 1ix=1,nx 
demi) ny 
Sum (1x7 1 0 d0 
enddo 
enddo 
‘ 
ele) ALS Laake 
Joe dt 1 1136 
Goma — ny 
P—cCaige loc )a 2) + Yi iy) * ¥2 ) 
ieremGee ten deaa=-20) Goto 110 


5,5 


110 


16.8) 


* 


* 


ern(1,1x,1ly)=.5d0*ar(n,27ix, ly) bowaee a ee 
(1.d0/rhol (1x, 1y) =1. 607 ene? aa ee 
ZO. (Thy, 1oe 1 y eizion  esiay a) 

ean(1,ix,iy) =.5d0*aa(n,1,1x, iy) hoa (ix/aea 
(1.d0/ rhe (a4) 1... d0/ wie ie ee 
ZO(m; ix, ay) * Zip (1 ace a 

erm(1,ix, iy) =.5d0*ar(m,1,1x,1y) 3) for (ike es 
(1.d0/rhol (1x,iy) -1.dG@/ rhoeZ (122731 ee 
yOu eabp.e tse) wevalor IL. those 147) 

eam(1,ix,ly)=.5d0*aa(m,1,ix,iy) + hpa(ix,iy)* 
(1.d0/rhol (1x, iy) -1.d0/ rho2{ 15a ae 
ZO (Mpa LV) * Zine, loca aie oe 

continue 

enddo 
enddo 
enddo 


do 1x=1,nx 
dou miy= ay, 
ler (hail, saben 
sum(ix,1ly)=Sum(ix,iy) + 
ern (15671) * erm (i oo ee 
ean(1,ix,iy) *eam(1,ix, iy) 
enddo 
enddo 
enddo 


if (n.eq.np.and.m.eq.mp) then 

write(6,*)"’ sum form im, in eae 

WrIte(6, LOO) s (Suma x, tec enh — 2 lee 
EOEMar(S ( loerel2 . 3) 
endif 


do ix=1,nx 
Ome aay, 
armn (15¢,4'y)=arim my tee eae) 
aamn(ix,iy)=aa(m,n,ix,iy) 
enddo 
enddo 


call partial (armn,nx,ny,x,y,arpr,arpa,nw, 
arpx,arpy,c, wk) 
call partial (aamn,nx,ny,x,y,aapr,aapa,nw, 
aapx, aapy,c, wk) 


Glow abe il, bt 
Gon dy— eng, 
r=dsqrt (x(1x)**2 + y(1y) **2) 
1£ (rc. lec ive 20) “ener 
bl (ix, 1y)=.5d0*arpr (130) 10 eee ee Ly) 
(1.d0/rhol (Gay iy) 4h G07 ietie7 alse yess Zion ney ae 


Se 


+ + + + 


end 
enddo 


return 
end 


Pinoignaia iy) Zlotor (ix, 1) ) 


else 

gw ge SG Owonapra(aac ie 5A0*ar (m,;nl,1x, ly) /r+ 
.5*d0*aapa(ix,iy)/r - sum(ix,iy) - 
(ileal Oy) Gods | toca) 
ieaOyTne2 (aoe, 1) ) ~Zoim.1x pty) * 
‘(Mey@he cual, SL us Al@hay ee | alo.Ge aliees 
Rpaidec aly = zbiioa (iss, iy) /r**2) 

endif 

do 


Sy 


wm Km KK KKK KKK KKK KKK KK KKK KKK KKK KKK KKK K KK KK KKK KKK KKK KKK KKK KKK KKK 


* * 
* This subroutine computes the partial derivatives with * 
* respect to range and azimuthal angle of a given function * 
* £(x,y). It uses a bicubic spline to calculate the * 
* cartesian derivatives and then perform a coordinate * 
* transformation to cylindrical coordinates. * 
* * 
KEKKKKKKKKKKKK KK KKK KKK KK KKK KKK KKK KK KKK KKK KKK KK KKK KKK KKK KKK KK 
Cc 
subroutine partial (£,nx,;ny,x,y, Epr, Epa; maweee see 
* CG vies 
Cc 
implicit real*8 (a-h,o-Z) 
real*8 £ (nx,ny),x(mxX),y (ny) 7 fpr (ne ny) fbea meee 2 
real*é@ wk(nw) , fpx (nx, ny) | tpy (ax, ny) Se (ee 2g 
Cc 
external ibcccu 
Cc 
lie =o 
pi=dacos(-1.d0) 
S 
call ibeeccut i -x7mx, yy Ce wie ore) 
e 
Gia e dsc lpnioc 
dor aly— ay, 
EO 1 sy, Cit 2 ps placings) 
Epy (2X21 a) =e ae een 
LE (x(1x)2eq.0-d0) srhen 
theta=dsign (y (iy) ,1.d0) *pi/2.d0 
else 
theta=datan (y (iy) /x (ix) ) 
endif 
P=GSqQrE ¢( 1) 842 437) ee) 
Epr (ix, ly) =fpx(ix,21y) “deos (enera) + 
* fpy (ix,1iy) *dsin(theta) 
fpa (ix, ly) =*Epx (ix, ly) *r*dsin (theta) + 
* Epy (ix, ly) *r*dcos (theta) 
enddo 
enddo 
S 
IgiSe (UliG ial 
end 
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